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Abstract We study the nearest neighbors one dimensional uniform q-model of force 
fluctuations in bead packs [Coppersmith et al (1996)], a stochastic model to simulate the 
stress of granular media in two dimensional silos. The vertical coordinate plays the role 
of time, and the horizontal coordinate the role of space. The process is a discrete time 
Markov process with state space R'^-'^'-'-'^j'. At each layer (time), the weight supported by 
each grain is a random variable of mean one (its own weight) plus the sum of random 
fractions of the weights supported by the nearest neighboring grains at the previous layer. 
The fraction of the weight given to the right neighbor of the successive layer is a uniform 
random variable in [0, 1] independent of everything. The remaining weight is given to the 
left neighbor. In the boundaries, a uniform fraction of the weight leans on the wall of the 
silo. This corresponds to absorbing boundary conditions. For this model we show that 
there exists a unique invariant measure. The mean weight at site i under the invariant 
measure is i{N + 1 — i); we prove that its variance is |(i(A + 1 — i))^ + O(A^) and the 
covariances between grains i ^ j are of order 0{N^). Moreover, as A — oo, the law 
under the invariant measure of the weights divided by A^ around site (integer part of) 
rA, r G (0, 1), converges to a product of gamma distributions with parameters 2 and 
2(r(l — r))~^ (sum of two exponentials of mean r(l — r)/2). Liu et al (1995) proved that 
for a silo with infinitely many weightless grains, any product of gamma distributions with 
parameters 2 and 2/p with p G [0, cxd) are invariant. Our result shows that as the silo 
grows, the model selects exactly one of these Gamma's at each macroscopic place. 
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1 Introduction 



In the nearest neighbors one dimensional uniform q-model of force fluctuations in bead 
packs the grains are arranged in layers inside a silo; each grain has random weight with 
mean 1 and supports part of the weight of the grains of higher layers. The weight sup- 
ported by each grain (plus its own weight) is randomly distributed between the two 
neighboring grains of the following lower layer. The fraction of the weight given to the 
right nearest neighbor is chosen uniformly in the interval [0, 1]; the left nearest neighbor 
receives the complementary weight. The walls of the silo absorb a random uniform frac- 
tion of the weight supported by the boundary grains. The resulting model is described 
by the following evolution equations. Fix a silo width and for z G {1, . . . , A^}: 

Wl'ii) = V^(z) + W,^i(z + l)?7i_i(2 + l) + W^,^i(z-l)(l-t/t_i(z-l)), (1.1) 
Wl^{Q) = Wt^{N + l) = 0, (1.2) 

where Vt{i) is the random weight of grain i of level t and Ut{i) is the random fraction of 
the weight grain i of layer t gives to grain i — 1 of layer t + 1. The boundary conditions 
( |1.2|) represent the absorbing wall. The random weights Vt{i) have mean one: MVt{i) = 1. 
The process W^^ is Markov in the state space (M"*")^. We are interested in stationary 
measures /x^ for this process. Under a stationary measure, the mean weights w'^{i) must 
satisfy the following equations. For i G {1, . . . , A^}: 

w'^it) = l + iw;^(z-l) + ^w7^(« + l), 

w;^(0) = w^'iN+l) = 0. (1.3) 

obtained by taking expectations in ( |1.1|) and assuming that the law at time t is the same 
as the law at time t — 1. These equations have a unique solution, that turns out to be 
quadratic: 

w^{i) = i{N + l-i). (1.4) 

This implies that all candidates to be stationary measures must have the same mean profile 
( |1.4j) . Our first result uses (|1.4|) to show the existence of an invariant measure. Indeed, 
starting with any measure with this quadratic profile, the means will be the same at all 
times. This guarantees that the sequence of measures indexed by times t G {0, 1,2,.. .} 
is tight {i.e. the mass does not escape to infinity; see the formal definition after display 
( pl2|) below) and by taking Cesaro limits along convergent subsequences, one can show 
that there exists a stationary measure. To prove uniqueness we realize simultaneously 
two versions of the process starting with different initial invariant distributions using the 
same sequence Ut for both evolutions. This is called coupling in the probabilistic literature. 
Under the coupling all weight added to the system after time zero evolves identically in 
both versions, while the discrepant weight present at time zero is distributed among the 
neighbors to eventually get lost at the boundaries. As a consequence, both versions will 
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have the same distribution in the hmit as t — > cxd. This is the key ingredient to show that 
there is a unique invariant measure fi^ for the silo. 

Next step is to study properties of the invariant measure /i^. It seems hard to exphcitly 
describe the invariant measure for each finite N, but it is possible to study the asymptotic 
behavior when grows. We show that the covariances between the weights supported by 
different grains under the invariant measure are bounded by a constant times and that 
the variance at site i differs from {w^ /2 by at most a constant times A^'^. This implies 
that if one normalizes the stationary weights (distributed according to /i^) dividing 
by A^^, one obtains, in the limit when N ^ oo, that the covariances of the normalized 
weights vanish and that the variance of {[rN])/N'^ converges to (r(l — r))^/2. Here 
[rN] is the integer part of rN. 

We also show that the law of the vector /N"^ around site [rN] converges to a 
product of gamma distributions with parameters 2 and 2/(r(l — r)), the sum of two 
independent exponentials of mean r(l — r)/2 each. To understand this result divide 
( |1 . 1| ) by iV^ and take N to infinity to get that the limiting values must satisfy the same 
equation but in the infinite lattice Z and without the Vt (weightless grains). The limiting 
equations correspond to the evolution of a model that we call infinite silo model (ISM). 
As Coppersmith et al (1996) did for the periodic computation shows that any 

space-homogeneous product of gamma distributions is invariant for the ISM. The key 
observation to show this is the fact that if X and Y are independent exponentials with 
mean p and U is independent of X, Y and uniformly distributed in [0, 1], then U{X + Y) 
and (1 — U){X + Y) are again independent exponentials with mean p. Actually we show 
more: following an approach of Liggett (1976) we prove that all invariant measures for 
the ISM that are also shift invariant are mixtures of products of gammas with parameters 
2 and p, with p e M+. Hence, the distribution of the vector (W^{[rN] + i)/N'^ : i E Z) 
converges to a mixture of products of gamma distributions, as A^ — > oo. Since the expected 
value of Vr'^([rA^] + tj/N"^ converges to r(l — r) for all £ G Z, it remains to identify its 
limiting distribution as the gamma with parameters 2 and 2/(r(l — r)). This is done 
using the fact that the normalized covariances vanish as A^ — > oo, as a consequence of our 
previous computation of the limiting covariances. 

Several authors proposed mathematical models for silos. The model in which a grain 
lies its weight on to the lower neighbors was introduced by Harr (1977) and explored by 
others, for example, Liu et al. (1995) introduced the model defined by ( ^.1[ ). The model 
with zero boundary condition was studied by Peralta-Fabi, Malaga and Rechtman (1997). 
Coppersmith et al. (1996) developed mean field computations for these models and in the 
uniform case, they conclude that the product of gamma distributions with parameters 2 
and 2/ p are invariant for a silo with periodic boundary conditions and weightless grains; 
they also obtained analogous results for dimensions d > 2. The main contribution of 
Theorem p.l2| is to prove that as the size A^ of the silo with zero boundary conditions 
increases, the invariant measure around a macroscopic site rN converges precisely to the 
gamma distribution with parameters 2 and 2/(r(l — r)). Rajesh and Majumdar (2000) 
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compute the space-space and space-time correlations for tiie periodic model starting with 
the zero configuration. Claudin, Bouchaud, Gates and Wittmer (1998) compute the 
spatiotemporal correlations in a continuum approximation when the mass of each grain is 
rcscalcd so that it goes to zero but Rajesh and Majumdar (2000) criticize the approach. 
Socolar (1998) shows numerical data for horizontal correlations for the model with periodic 
boundary conditions. Mueth, Jaeger and Nagel (1998) show experimental data for force 
correlation in a silo filled with glass beads. Essentially they conclude, in agreement with 
the q-model, that the correlations vanish. Lewandowska, Mathur and Yu (2000) computed 
the variances of stresses in the q model with weightless grains. Krug and Garcia (2000) 
and Rajesh and Majumdar (2000) study related conserved-mass models with periodic 
boundary conditions and show that in some cases the process has a stationary measure 
with non zero correlations. 



2 Definitions and results 

To define the model let {Ut) '■= {Ut(i) : i G Z, t > 0} be a family of independent uniform 
random variables in [0, 1] and (Vt) := {Vt{i) : i G Z, t > 0} be a family of iid positive 
random variables with mean 1 and variance a = YVt{i) < oo. Furthermore assume (Vt) 
and {Ut) to be independent families. 
Fix N > 1, consider the finite box 

:= {1,...,N} 

and denote Wf^{i) the weight carried by a grain located at the ith position at level t. Fix 
an initial configuration G [0, oo)^ and define inductively 

Wfi{} = Vt{i) + Wt''_,{i + l)Ut_,{i + l) + Wf_,{i-l){l-Ut-i{i-l)), iGA^ 
Wf{0) = Wt^{N + l) = 0. (2.1) 

Let Wl^ = {W^^ii) : i G A^); then {W^^ : t > 1) is a discrete time Markov chain on 
[0,oo)^ . Each grain j of layer t gives a fraction chosen uniformly in [0, 1] of its own 
weight plus the total weight it supports from the previous layers to grain j — 1 of the 
successive layer t + 1 (which we can think is below t) and the remaining to grain j + 1. 
The weight distributed to grains outside A^ is thought of as being absorbed by the walls 
of the silo at sites and -|- 1. 

For a measure u on [0, oo)^ let i'S'^{t) be the measure defined by 

uS'^it)! = EJ{Wf) := J u{dW)E{f{Wf)\W,'^ = W). 

where E and P are the expectation and probability defined with respect to the probability 
space induced by {Ut : t > 0) and (Vt : t >0). 
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We say that a measure //^ is invariant for the process if fi'^ S^{t) = /i^. 

If W^^ has an invariant measure /x^ its mean heights w^{i) := J fi^{dW)W{i) have 
to satisfy the following system of equations (taking expectations in (|2.1| ) and assuming 
that the distribution of does not depend on t): 

w'^it) = l + ^w^{t-l) + ^w''{t + l), for^GA^; 

w^{i) = 0, for z G Z \ . (2.2) 

Notice that the only finite solution w^{-) of ( p. 21) has quadratic profile: 

w^{i) = i{N + l-i), hi i = 0,...,N + 1. (2.3) 

(It is the expected time to exit A^ for a symmetric nearest neighbors random walk starting 
at i.) 

Our first result says that for each N the silo with absorbing boundary conditions 
admits a unique invariant measure and that the process starting with any initial condition 
converges to the invariant measure. 

Theorem 2.4 There exists a unique invariant measure fi^ for the process defined 
by ( ^■^ ). The mean weights under fi^ are finite and satisfy l\2.3j . Furthermore, for any 
initial configuration , the law of Wf^ converges to fi^ . 

Denote the covariances of weights i and j under fi^ by 

a^(z,j) := [ fi''{dW)W{z)WU) - w''{z)w''U). 



We show the following bounds for the covariances. 

Theorem 2.5 The covariances of fi^ are of order : there exists a positive constant C 
such that for all N 

\a^{i,3)\ < CN\ 1^3, t^jeAN (2.6) 



and 



a [1,1] 



< CN'\ igAn. (2.7) 



2 

Numerical analysis performed in Section 2 suggests 

a^([xAr],[yAr]) ^ _r_{x^ _ ^.^ a^([xiV], [xAT]) - ^ r(x,x) 

(2.8) 
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for x,y & [0, 1], X 7^ y, where r{x, y) is the unique weak solution of the PDE 

- An-a(a;,y)(^AM + 2^j = f{x,y) in = (0, 1) x (0, 1) 

u = in (2.9) 

where a{x, x) = 1/2 and a{x, y) = for x ^ y and A stands for the Laplacian operator. 
The right-hand-side is given by /(x, y) = 3\/2(6(x(l — x) — l)S{x — y), where 5 is Dirac's 
distribution. An approximate plot of the function r obtained numerically can be seen in 
Figure 1. The numerical solution also suggests that 

hm — ^ = c{t,j) > 0. (2.10) 



The meaning of ( |2.10|) for i = j = 1 is that the variance of the weight discharged to the 
left wall at each level is of order A^^, the square of the mean. Hence the fluctuations are 
of order A^, the order of the mean. By symmetry the same happens to the right wall. On 
the other hand, ( |2.10| ) for i ^ j means that the stationary correlations divided by the 
square of the mean stay strictly positive around the origin; this discards the hypothesis 
that W'^{-)/N'^ converges to a product measure. In contrast, the correlations divided by 
the square of the mean vanish in the interior of the silo as we show in the next theorem. 
For p > define 



product measure on [0, oo)^ with marginals Gamma (2,2/p) (2-11) 



(the sum of two independent exponentials with mean p/2 each). Our main result says 
that if W is distributed according to the invariant measure /i^, then as N —>■ oo, the law 
of W/N"^ around the site [rN] converges to Urii-r)- Let be the translation operator by 
z on Z and 9^ be the operator "divide by fc": for W e [0, cxd)^, {TiW){j) = W{j — i) and 

{QkW){j) = w{j)/k. 

Theorem 2.12 Let jj.^ be the unique invariant measure for the process W/^ . Then for 
each r G (0, 1) 

lim r[r AT] 9 Ar2 p^ = Uril-r) (2-13) 

weakly, where [■] is the integer part. 



The proof of Theorem p.l2| lies on two results. The first one is related with the limit 
as A^ — s> oo of the weights of the silo model divided by A^^ in sites around a macroscopic 
position rN. Its distribution converges to the law of another process in (M"*")^ whose 
evolution is the same described by ( p.l|) but with Vt{i) = — no extra mass is added in 
this system. This process is called the infinite silo model, ISM. We follow an approach of 
Liggett (1976) to show that in the nearest neighbors uniform case, the invariant and trans- 
lation invariant measures for the ISM are mixtures of product of Gamma distributions 
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(Proposition 4.11| below). Combining this with the fact that the stationary space-space 
correlations go to zero, proved in Theorem we can single out the product measure as 
the limit in ( |2.13|) . 

Theorems and |2.5| are proven in Section ^ A numerical study of the corrections for 
the correlations in Theorem |2.5| is given at the end of Section |^. Theorem |2.12| is proven 
in Section 4. In Section 3 we study the infinite silo model. 



3 The finite stationary silo 

In this section we study the silo for fixed width A^. We prove the existence of an invariant 
measure and study its correlations. 

Notice that ( |2.1| ) describes the evolution of two independent systems: {W^{i) : i + t 
even} and {Wf^{i) : i + t odd}. It is convenient to keep simultaneously track of both 
systems. 

Proof of Theorem \Z.4\ If we consider any initial distribution uq for , not necessarily 
invariant, satisfying Ej,q [1^5^(2)] = i{N + 1 — i), x G A^, then for alH > 

Ko[Wf{i)] = t{N+l-z),zeA''. (3.1) 

Since by Markov inequality 

P.J<(.)>M]< ^-°[^'^^^^ (3.2) 

the sequence {Wt : t > 0) is tight [this means that for all £ > there exists a A' > 
such that ¥{Wt(i) > K) < 1 — e for all i,t\. Tightness implies that the Cesaro means 
ut '■= {^/T)J2t=i^'^i^)^o converge through subsequences to an invariant measure (see 
Liggett (1976), for instance) called /i'^. 

If (tfc) is a subsequence such that ~^ 1^^ 1 as /c — 00, then there exist with law z/j^. 
and W with law /i^, all defined in the same probability space, such that for i = 1, . . . , A^, 
Wk{i) — > W{i) almost surely (Theorem 8.1, page 25 in Thorisson (2000); this is usually 
called "almost sure Skorohod representation for distribution convergence"). Since the 
means are independent of t, Fatou's Lemma (E(liminffc Wk{i)) < liminffcEiyfc(i)) implies 
that ¥W{i) < i{N + 1 — i), that is, the means under /x^ are bounded by the quadratic 
profile (p^). But if the means are bounded and satisfy the equations (p.2|), then /x^ must 
have the quadratic profile (|2.3|). 

In order to prove uniqueness we use coupling. Construct two processes Xt and 
satisfying ( p.l|) using the same uniform random variables {Ut{i) : i,t) and {Vt{i) : i,t). 
Assume the initial configurations satisfy 

Xo(i) > lo(«) for i e A^. (3.3) 
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Then it is easy to see that at any time t > 

Xt{i) > Yt{i) for i G (3.4) 
(in this case we say that the process is attractive) . Denote 

1 ^ 



N 
1=1 



We have 



Dit) = D{t - 1) - ^^i^[x,_i(i) - r,_i(i)] - ^^f^[x,_i(iv) - r,_i(iv)] . (3.5) 

The numbers D(l), Z?(2), . . . form a bounded below by non-increasing sequence of non- 
negative random variables. Hence D{t) converges almost surely. This implies that its 
increments must go almost surely to zero. This and ( |3.5| ) imply that 



[Xi_i(l)-Fi_i(l)]^0 and [X,_i(iV) - r,_i(iV)] ^ , 
almost surely, as t ^ oo. Looking at the equations 

Xt-i{l) - Yt-i{l) = (1 - Ut^2{2)) [Xi_2(2) - 1^-2(2)] 

and 

Xi_i(iV) - F,_i(iV) = Ut-2{N - 1) [Xi_2(iV - 1) - Yt_2{N - 1)] , 
it is immediate to see that 

Xi_2(2) - y;-2(2) ^ and X^.^iN - 1) - Y^^^iN - 1) , ^ 

almost surely, as t — oo. Using the same argument we have that for all i G 

Xtii) — Ytii) —^0, as t ^ oo , 

almost surely as t — oo. 

Let (^o(0 • ^ ^ ^^'i (W^o(0 • ^ ^ arbitrary random vectors (not necessarily 
ordered). Define 

Xo(z) = Zo{t) V Wo{i) and yo(^) = ^o(^) A Wo(0 • 

We have 

Y^{i)<W^{i)<X^{i). 

Couple the four processes using equation (|2.1| ). This means that the four evolutions 
are performed using the same random variables Vt{i) and Ut{i)- The previous argument 
implies 

WAi) - YA^ ^ and XAi) - WAt) ^ , 



8 



almost surely as t — > oo. Since 

Wt{i) - Yt{i) < \Wt{i) - Zt{i)\ < Xt{i) - Wt{i), 

we have 

Wt{i)-Zt{i)^Q, (3.6) 

almost surely as t — > oo. 

Hence, if we pick Wq from the distribution /i^ and Zq from an arbitrary invariant 
distribution z/, we have This proves that /x^ is the unique invariant measure. 

The limit ( |3.6| ) also shows that the process starting from an arbitrary configuration 
converges to /i^. □ 



Proof of Theorem |2.5| At the end of the proof we show that < oo for all i, j 

and A^. If is distributed according to the invariant measure, so is , and one can 
use (|2.1| ) to show that a^{i,j) satisfies the system of equations: 

O O t: 

+ \^''{^ + 1,^-1) + ^iw'^i^ + l)r + ^iw'^i^ - l))^ ^ G 
a^{i, i + 2) = ^a^{i + l,i + 3) + ^a^{i - 1, i + 3) + ^a^{i - 1, z + 1) 

+ V(z + 1,^ + 1) - + ^e{l,...,N-2} 

a^{t,t-2) = + + 1,^-3) + + 1,^-3) 

+ 1,^-1) - ^e{3,...,Ar} 

+ V(z + l,j-l), z,jeA^, K-j|>2 

= 0, zG{0,Ar+l}orjG{0,Ar+l}. (3.7) 
Define, for i e {0, 1, . . . , A^ + 1} : 

Rii,i) = -2a^it,t) + {w''ii)f, 

R{i,3) = -3(7^(2, j) for 2^ J. (3.8) 
The system ( |3.7|) translates into the following system for the matrix R: 

Riz,j) = Ki^)l{^=J}+J2piihJ),i^,f))Ri^,f) ^,jeA^, 

R{i,j) = i e {0,Ar + l} or J G {0,A^+1}, (3.9) 
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where 



= 6i{N + - (N + lf - {l + 2a) , (3.10) 



for i G A and 



1/4, if i 7^ j and |i — i'l-lj — j'l = 1 



otherwise. 



for i,j G {1, . . . , A^} and 



I- otherwise 

if i G {0, 1} or j G {0, 1} (this means absorbing borders). Let Xn'"*^ be a Markov 
chain with transition probabihties p{-, ■) and initial state Defining 

N 

rp{i,j) ■=J2Y1 l{^n '^'^ = ^)} (3-12) 
n>0 ^=1 

we have that ET^*'-^^ is the unique finite solution of (|3.9|) . Assuming that \a^{i,j) \ is finite, 
then \R{i,j)\ is also finite and R{i,j) = ET'^^'^l 

To get an upperbound for R{i,j) we consider a random walk X^*'"*^ on the diamond 

:= {(z,j) G : i + J even, 2 < i + j < 2X, |i - j| < X + 1} 

with transition probabilities 

Piihj), («', j')) = Piihj), {^'d')) given by (3.11) 

if (i,j) G {(2,j) G A^ X A^ : 4 < i+j < 2X-4, |i - j| < X - 1} (3.13) 

(the interior of the diamond); periodic boundary conditions along two of the sides of the 
diamond: 

p((z,j),(A^ + l-J,Ar + l-i)) = p((iV + l-j,iV + l-z),(z,j)) = 1/4 

if i+j = 2, i^i 

p{{^,J),{^',f)) = p{{N + l-j,N + l-z),{N + l-f,N + l-z')) = 1/4 
if t+j=2, t^j, |^-^'|. 1^-/1 = 1, («',/)gA^ 
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p{{N,N),it,j)) 



1/3 if {(2,2), (iV,iV)} 

1/6 if(i,j)e{(2,0),(0,2)} 

ri/3 if(i,i)e{(l,l),(7V-l,7V-l)} 

\ 1/6 if (i, i) e {(TV - 1, TV + 1), (TV + 1, TV - 1)} 



and absorbing probabilities along the other two sides: 

Pii^,J),i^,J)) = 1, (^,j)eA^, |z-j|=TV + l. 



(3.14) 
(3.15) 

(3.16) 



The point of this random walk is that it is "one dimensional" in the sense that the 
one-dimensional process defined by Z'-^ := (xt^^)i - {xt% (the difference of the 
coordinates of Xn'^^) is Markovian and has transition probabilities 



Po{i,j) 



'1/4 ifi ^ {-Ar-l,0,iV + l} and = 2 

1/2 if J =i and i ^ {-A^- 1,0, A^+1} 

2/3 if i = j = 

1/6 if i = and and |i — j | =2 

1 if i = j and i e {-A^- l,iV + 1} 



(3.17) 



for even i e {-N - 1, . . . , iV + 1}. 

Define K = max^ K{£), K = min^ K{£). We have I{ < < K and both K, K are of 



order of N . Define 



ih3) 



N 



r 



n>0 l=\ 

N 



n>0 e=i 

A simple coupling argument shows that for all A^, 

rp(i,i) ^ J^iiJ) ^ <; Y'^*'-?) ^ Tj^i-hi) 



(3.18) 



On the other hand T^*'*^ and J"^*'*^ do not depend on i and can be expressed in function 
of a geometric random sum of geometric random variables: 



(3.19) 



fe=i 



fe=i 



where ¥{S > n) = (^)""S n > 1 and P(Mfe > m) = (2/3)"'-\ m > 1. Furthermore 
(Mfc) arc iid and independent of S. The randoni variable S counts the number of times 
the chain Z„ starting at 2 or —2 comes back to zero before arriving toA'^ + lor— A?^ — 1 
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while Mk is the time the chain Z„ spends in its kth. visit to the origin. Taking expectations 
in ( p.l9| ) and using Wald identity we get 

e¥'''^ = -{N + 1)K ; ET^'''^ = -{N + 1)K. (3.20) 

Taking expectations in ( |3.18| ) and using the fact that both ^ and K are of order A^^, we 
have proved that if R{i,j) is finite, then 

- CN^ < R{i,3) < R{i,3) < R{i,]) < CN^ . (3.21) 

This, (^) and (^ show (^ and 

Finiteness of a^{i,j) Let be a non negative random configuration chosen accord- 
ing to a product measure with marginals satisfying 

YWo^'it) = {w^{t))V2 (3.22) 
(the covariances are null). The covariances at time t are given by 

^a^,J)- i _i^0(,^^-) ifz^j' ^'^■^'^^ 



where R^{i,j) evolves according to 



R',it,j) = l{^ = + J2 Piihj)^ {i\3'))RlS'.3) , (3.24) 

with i?Q(z, j) = for i, j G {0, . . . , + 1}. Let Rt and be the systems defined in the 
diamond with initial condition RQ{i,j) = R^{i,j) = 0, i,j G and evolving with 
the equations 

:R°(^,j) = l{^=J}K+J2PiihJ),{^,f))'Ift^^{^',f)■ (3.25) 
^S.J) = ■^{^=j}K+Y,Piihj)S'.j'Wt-S'.j')- (3.26) 



It is clear that 



On the other hand 



RAhj) < Rtihj) < CN\ (3.27) 
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where Rt = R is the stationary solution of (|3.25|) . That is, R^{i,j) is bounded by the 
stationary solution of (|3.25|) which in turn is bounded as in ( p.21| ). The same argument 
shows that \R^{i,j)\ is uniformly bounded in t. Hence \R^{i,j)\ is uniformly bounded 

we need to take 

,N 



To show the finiteness of a 



Ni 



in t and by (|3.23|) so is \(Tf{i,j] 

the limit as t — > oo. Since the distribution of Wf^li), converges weakly to /i^, as 
t oo, we can again use Theorem 8.1, page 25 in Thorisson (2000) to define Wl^{i) and 
ly^(z) in the same space with marginal laws z//^ and yU^, respectively, such that 



tN I 



converges almost surely to W 



N 1 



as t ^ oo. Now, a 



N , 



1,1) 



NW 



N I 



{mV^{i)f. Since ¥W 



N I 



N I 



E(liminft(W,^(i))2) < liminftE(iy^(^))2 



for all t, Fatou's Lemma (which in this case says 
implies 



< liminf ai^fi, z) < 



[w 



by ( p.23|) and (|3.21|) . This settles the finiteness of the diagonal covariances. By 
the other covariances must be finite. □ 



(3.28) 

), also 



'Solution_N=127' + 



0.3 p 
0.25 - 




Figure 1: Graphic of the computed solution r = R/ {N + 1)^ for = 127. 



Numerical Analysis of the equation for the correction R We carried out a numer- 
ical investigation about the solution of equation (|3T9|). This equation can be seen as a dis- 
cretized version of the partial differential equation ( |2.9|) in the following way. We first di- 
vide system ( p.9| ) by (A^+1)^, obtaining a system for the new variable r^j := Rij/ (A^+1)^, 
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0.3 




-0.1 I 1 1 1 1 1 

0.2 0.4 0.6 0.8 1 



Figure 2: Graphic of the computed solution r = R/{N + 1)^ over the diagonal (z = j) 
for several values of N. 



which is to be seen as an approximation to the value of a function r{x,y) defined on 
Q = [0,1] X [0,1] at the point {xi,yj) = {ih,jh), where h = 1/(A^+ 1)^. The equations 
for i j are of the form 

4rjj — Tj-ij-i — Ti-ij+i — rj+ij_i — ri^ij^i = 

which is a discrete form of Ar = 0, where A denotes the (rotation invariant) Laplace 
operator. Now, for i = j we multiply ( p.9| ) by 3h = 3/{h'^{N + 1)^) and obtain the 
equations: 

1 1 

= 3^2 (6x,(l - x,) - 1 - (1 + 2a)/i2) - i) , (3.29) 

where = ih and is a piecewise linear approximation to Dirac's distribution, as- 
suming the value l/(\/2/i) at the origin, and being zero at the other points of an equally 
spaced mesh (with mesh-space H = \/2h). Now, we observe that the first term in equa- 
tion ( |3.29| ) is a discretization of minus the Laplacian of r (in rotated form) and that 
the second expression corresponds to minus the second derivative of r along the diagonal 
X = y (which can be written as — l/2Ar — d'^r /dxdy). The right-hand-side of ( p.29| ) 
approximates the function 3y/2{6x{l — x) — l)5{x — y) over the diagonal. All together we 
obtain equation ( |2.9| ), which shall be interpreted in the weak sense. 
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To compute the solution of equation ( p.9|) we have developed a multigrid method 
composed by red-black Gauss-Seidel relaxations, full- weighting of residuals for fine-to- 
coarse mesh transfers and of bilinear interpolation for coarse-to-fine corrections (see eg. 
Hackbusch (1985) for definitions). The scheme is very fast, it takes only a few seconds 
to solve equation ( p.9|) for = 511 on a Sun- workstation and around 15 minutes for the 
solution for = 4095 on a Dec-alpha (the equation in this case has more than 16 million 
unknowns). The equations were solved to the level of machine precision, using double 
precision (16-byte) words. 




Figure 3: Graphic of the values of -R(l, 1)/(A^ + 1)^ for several values of log2(A^ + 1). 



The results shown in Figures 1,2,3 were obtained with a = 0, but the same asymptotics 
are expected for other values of a. In Figure 1 we present the function r = R/{N + 1)^, 
with a resolution of = 127 (the function doesn't change much visually for larger values 
of N). In Figure 2 we plot r = R/{N + 1)^ over the diagonal {i = j}, for several values 
of up to 511 - for larger values the graphics are indistinguishable. We can observe 
the convergence of r as A^ grows (to the weak solution of the partial differential equation 
( |2.9| )). The upper and lower constant bounds for R/ (A^-|-l)'^ we proved to exist in theorem 
( p.5|) , are attained over the diagonal {i = j} and their approximate values can be read 
from Figure 2. 

Another aspect we tried to access from the numerical solutions is the behaviour of 
R{1, 1)/(A^-|- 1)^ (which corresponds to the derivative of r over the diagonal - the function 
of Figure 2 - at x = 0). We computed these values for A^ = 15 up to A^ = 4095 (with 
A^ 4- 1 as powers of two). The results are presented graphically in Figure 3, where A^ is 
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displayed in log scale. The results suggest that as — > oo 1)/(A^ + 1)^ 
a constant (around —3). This together with the definition of given in (|3 
( |2.10|) . The constant c(l, 1) seems to be around 4. 



converges to 
supports 



4 The infinite silo model 



To prove Theorem p.l2| we will use a spatially infinite limit of the silo model. This model 
represents the transmission of weight among weightless grains in an infinitely large silo. 
For each i G Z let ?7t(i) denote the weight supported by site i at time t > 0; the dynamics 
is represented by the system of equations: for i G Z and t>l such that it is even, 

r/i(z) = r/i_i(z - l)Ut^i{i - 1) + r/i_i(^ + 1)[1 - Ut-i{i + 1)] , (4.1) 

where {Ut{i) : ? G Z, t > 1) are i.i.d. uniform random variables in [0,1]. We call this 
process infinite silo model (ISM). Let Up be a product measure on [0, oo)^ with marginals 
Gamma(2, l/(2p)); p is the expected value of each marginal. Coppersmith et al. proved 
that the Up is invariant for the process with periodic boundary conditions. We show now 
that Up is reversible (and hence invariant) for the infinite system. 

Proposition 4.2 For each p > the measure Vp is reversible for the ISM, that is, for all 
cylinder functions f and g, Kpifivo)9iVt)) = ^•ypi9iVo)fiVt))- 

The proof is based on the following elementary lemma. 

Lemma 4.3 Let W be a random variable with distribution Gamma (2,p/2). Let U be 
a uniform random variable in [0, 1] independent of W . Then, UW and (1 — U)W are 
independent exponentially distributed random variables with mean p/2. 

Proof of Proposition |4.2| It suffices to prove that if tjq has law Up then, for cylinder 
functions / and g, 

^f{vo)givi) = %(^o)/(^i) . 

Let 

X(z) = Uoit)vo{^ ; Yit) = (1 - Uo{t))vo{^ . 

Under Up, {rio{i) : z G Z) is a family of iid gammas random variables. Then, by the 
lemma, {X{i) : i G Z) and {Y{i) : i G Z) are independent families of iid exponential 
random variables of mean p/2 and 

^f{Vo)giVi) = E/(X + Y)giX' + Y') , (4.4) 

where {X'{i),Y'{i)) := H{X,Y){i) := {Y{i - l),X{i + I)). Since H{H{X,Y)) = {X,Y) 
and {X',Y') has the same law as {X,Y) (iid exponential random variables), (f4.4| ) equals 
to 

EfiX' + Y')g{X + Y)= Eg{r]o)f{vi) ■ □ 
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One of the tools to prove Theorem 2.12 is to show that all invariant and translation 
invariant measures for the ISM can be written as convex combinations of those product 
measures. 

For any two initial configurations rj and ^, we couple the two processes rjt and using 
the same uniform random variables to update the process. Notice that this coupling is 
attractive, that is, 

V{i) < ^(0 for all i G Z, then r]t{'i} < ^t(z), for alH G Z. (4.5) 

Through out this section X and T (X and T) will denote, respectively, the set of invariant 
and translation invariant measures for the process (coupled process) defined by ( |4.1| ). 

Proposition 4.6 If u El (iT, then 

H{V,0;V>^ orr]<^} = l. (4.7) 

The proof of this proposition is based on the following lemma. 

Lemma 4.8 Let X , Y be identically distributed real random variables defined in the same 
probability space. Let U and V be identically distributed and independent random variables 
in [0,1] with E?7 = 1/2; U and V are assumed to be independent of X and Y. If Z = 
UX + VY IS distributed as X (andY), then P(XF > 0) = 1. 

Proof. From the definition of Z, 

E(Z1{Z > 0}) = ^{UX1{Z > 0}) + E{VY1{Z > 0}) 
= 2E{UX1{Z > 0}) 

= 2E{UX1{Z > 0, X > 0}) + 2E{UX1{Z > 0, X < 0}) 
= 2E{UX1{X > 0}) - 2E{UX1{X > 0, Z < 0})}) 
+ 2E{UX1{Z > 0,X < 0}) . 

The second identity follows from the fact that UX1{Z > 0} and VY1{Z > 0} have the 
same law. The third and fourth identities come from set operations. Since X and Z are 
identically distributed and U is independent of X, E{Z1{Z > 0}) = 2E(f/Xl{X > 0}). 
Canceling those terms we get 

= - E{UX1{X >0,Z < 0})}) + E{UX1{Z > 0, X < 0}) . 

This implies 

P(X>0,Z<0) = P(X<0,Z>0) = 

p(r>o,z<o) = p(r<o,z>o) = o (4.9) 

because (X, Z) and {Y, Z) have the same distribution. Hence 

p(x > 0, F < 0) = p(x > 0, r < 0, z > 0) + p(x > 0, r < 0, z < 0) = . 

The same argument shows that P(X < 0, F > 0) = 0. □ 
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Proof of Proposition 4^ Let (?7o;^o) be distributed according to z/ G X n T. Then 
the variables 

X = - 1) - r7o(2 - 1); Y = + 1) - r^,{t + 1); Z = U^) - Vii^) (4-10) 
satisfy the hypothesis of Lemma [4.8| . Hence, for all i G Z, 

1/(^(2 - 1) = v{i - 1)) = 1 or z/([eo(2 - 1) - - 1)] [eo(« + 1) - + 1)] > 0) = 1 , 
which implies the proposition. □ 
Proposition 4.11 If fi e {I HT) then 

fi = j Upd\{p), (4.12) 

where X is a probability measure on (0, 1). 

Proof. For each p > we can construct Up & I CiT, a measure on [0, oo)^ x [0, oo)^ with 
marginals /i and z/p (Liggett (1976)). 

We claim that /i G X fl T concentrate its mass on configurations with asymptotic 
density. That is, for a configuration 77 define 

1 

Mniv) = E ^(*) • (4.13) 

i=—n 

Then, for all a G [0, 1], 

H{r];\immi Mn{r]) < a< limsup M„(?7)} = 0. (4.14) 

In fact, if ( |4.14D is false, there exists an a such that 



i^aiiv^O'^V >^ or T] < > liminf M„(r7) < a < limsup M„(?7)} > 0, (4.15) 



which contradicts Proposition Therefore, taking = po < Pi < • • • < Pm = 1, with 
Pi = i/m, with m positive integer, we have 

m 

fI{■)=J2f^i■\C^)p^, (4.16) 

i=l 

where Q := lim„ Mn{r]) G {pi-i, Pi]} and pi = p{Ci), for i = 1, . . . , m. 
Define 

m m 

P„(-) = Xl^"-!^^' J^mi-) =^l'p^P^■ (4-17) 

1=1 i=l 
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By definitions (|4.16|) and ( [4 .171) and Proposition ( |4.6|) we can couple /x, fi and and show 



!—m 



for all m, where <st is the usual stochastic domination of measures defined by: fi <st if 
and only if for all non decreasing function /, fif < fi'f; equivalent ly, <st fJ'' if and only 
if there exists a joint realization {W, W) with marginals fi and fi' such that W < W. 
Defining A(p) = fi{ri : lim„M„(?7) < p) we have 



lim fi = lim jj,^ 



This proves the proposition. 



□ 



5 Proof of Theorem \2.12 



In this section we use the characterization of the set of invariant and translation invariant 



measures for the ISM given by Proposition [4.11| and the bounds on the covariances ( |2.6| ) to 
prove Theorem |2.12| . Let be the unique invariant measure for the silo in the finite box 
A^. We extend the definition of the measure /i^ in [0, oo)^ by setting fi^{ri : ri{i) = 0) = 1 
for all i^A^. 

Let Ut{i) be the uniform random variables defined in the introduction and define Ut 
as the matrix with entries 




(5.1) 



Define the truncations 



[Vfi^) 



V{i)l{i G A^}. 



(5.2) 



The evolution equation ( |2.1| ) is then a product of a left vector and a matrix plus an 
independent vector: 



Notice that the definition of [U^] guarantees W^f (j) = for j ^ A 
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Existence of the limit Let W'^ be a configuration distributed according to tlie invari- 
The maximum of ( p.3[ ) is attained for i G {N/2, [N + l)/2}, hence 



ant measure 



E 



iy^([rA^] +£) 
iV2 



< 



O 



2NJ 



(5.3) 



for all i G Zi. Therefore, the sequence of probability measures {T[rN]QN^I^^,N > 1} (that 

is, {law of ^ ^j^f^^"^'"* , N > 1}) is tight. Consider a convergent subsequence {Nk,k > 1} 
and let /i be the weak limit of this subsequence: 



II := lim T[rNk]^N^^ 



Nk 



(5.4) 



Invariance of the limit The equation /x^ = S^{t)iJi^ is equivalent to 

E/(W^^) = E/(iy^[U]^ + [Vf) , (5.5) 

for cylinder bounded /, where U has the same law as U^, is a vector with the same 
law as Vj+i and , U and V are independent. 

If 7] is distributed according to z/, then the invariance of v for the infinite silo model is 
equivalent to 

E/(r/) = E/(r/U) , (5.6) 

for cylinder bounded /, where rj and U are independent. 

Let fj^ = T[rNk]'S)N^W^'', = T[rNk]V and U'' = r[^Ar^]U, where {Ti'U){iJ) = \J{i - 
i,j — £). Take k so large such that the support of ^rN^lf is inside A^*^ \ {1, N^} (and the 
truncations are unnecessary). Then the law of fj^ must satisfy 

Efif) = Efiflj'' + (VViVfc^)) . (5.7) 

By translation invariance, U'^ and U have the same law, and V have the same law 
and by ( |5.4| ) f/^ converges in distribution to /i as — oo; hence /i must satisfy (^]6|), that 
is, e X. 

Translation invariance of the limit Let {W^'^, W^^ , W^^^) be identically null and 

iWiii^Wli,,Wlir) (5.8) 

for t > 0. Then Wf-^ < W/^ < Wf+^ for all t. By Theorem U (W^f "^ Wf, Wl^+^) 
converges in distribution to {W^~^, W'^ , W^~^^) a vector with marginals ri/x^~^, Tifi^ 
and also satisfying W^'"^ <st <st W'^^'^, where X <st Y means that the law 

of X is stochastically dominated by the law of Y (see the definition at the end of Section 
4). This implies 

. (5.9) 
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Analogously, the coupling 

{wii-,\w,%w,ir) (5.10) 



t+l 5 ^^t+l^ ^^t+l 



shows that 

ri/i^-' /i^ (5.11) 
From the coupling ( p.lO| ) and by Chebyshev, 

/ ^^+2(,) -iyiV(,) X ^ 2(Ar + 3-z)-z(Ar + l-^) ^ 2z 

ly^(i) - - 1) \ i(^N +l-i)-{i-l){N -i) _ N 



as — cx) for all i G {0, 1, . . .}. Therefore, if ^^''TrN^Qj^^2 converges weakly to ^ the 
same is true for fi^''~'^TrN,.+iQ f^^'^'^^'^rNk^N,,^- This and ( |5.9| ) imply ri/i = /i and 
hence that G T. 



Correlations of the limit. Since G X fl T, by Proposition |4.11| /i is a convex com- 
bination of gamma distributions with mean r(l — r). The covariances of n equal the 
variance of the density of fi: for all i j, 




^{drj)r]{i)ri{j) - ( / fi{dr])rj{€) ] = j J Up{dr])ri{i)r]{j)X{dp) - i^J pX{dp) ^ 

p'^Xidp) - (^J p\{dp)y > 0. (5.12) 

Since the joint distribution of the couple (^^^^^^^, '^^'^^'^t''^ ) (under p^'') converges 
in distribution to the joint distribution Q of {j]{i),rj{i)) (under p), then, as in the dis- 
cussion following ( ^.2| ), by Fatou's Lemma and Skorohod representation for distribution 
convergence, / xydQ < liminffc^oo/ xy Qk{x,y). Since by ( ^.3|) the means converge and 
limfc^oo " t^i''^ = (by Theorem |2.5|) the covariances of p are nonpositive. 



Conclusion: We have proved that p is a convex combination of Up with mean r(l — 
r) and with zero correlations. This means that p = Vr{i-r)- Since each convergent 
subsequence of {p^ TrisiQ^'^i A^ > 1} converges to Vr(i-r), so does the sequence. □ 
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